module utils

    use, intrinsic :: iso_c_binding
    use basic
    use response
    use spectra
    implicit none

contains

    subroutine r(acc,n,dt,zeta,P,ra,rv,rd,SM) bind(c)
        integer, intent(in) :: n, SM
        real(8), intent(in) :: dt, zeta
        real(8), intent(in) :: acc(n), P
        real(8), intent(out) :: ra(n),rv(n),rd(n)

        select case (SM)
            case (1)
                call rfreq(acc,n,dt,zeta,P,ra,rv,rd)
            case (2)
                call rnmk(acc,n,dt,zeta,P,ra,rv,rd)
            case (3)
                call rmixed(acc,n,dt,zeta,P,ra,rv,rd)
            case default
                call rmixed(acc,n,dt,zeta,P,ra,rv,rd)
        end select
    end subroutine r

    subroutine rnl(acc,n,dt,zeta,P,ra,rv,rd,SM,mu) bind(c)
        integer, intent(in) :: n, SM
        real(8), intent(in) :: dt, zeta, mu
        real(8), intent(in) :: acc(n), P
        real(8), intent(out) :: ra(n),rv(n),rd(n)

        integer :: i, model, maxiter
        integer, parameter :: npara = 8
        real(8) :: w, k0, fy, uy, rk, k1, muc
        real(8) :: pkd, pkf, para(npara), tol

        model = 0
        maxiter = 1000
        tol = 1.D-9

        rk = 0.D0
        w = pi2/P
        k0 = w*w
        k1 = k0*rk

        call r(acc,n,dt,zeta,P,ra,rv,rd,SM)
        pkd = abs(peak(rd,n))
        uy = pkd/mu
        pkf = pkd*k0

        if ( SM == 0 ) then
            para = [ k0,k0*uy,k1,0.D0,0.D0,0.D0,0.D0,0.D0 ]
            do i = 1, 1000, 1
                call rnmknl(acc,n,dt,zeta,P,ra,rv,rd,para,npara,model,tol,maxiter)
                pkd = abs(peak(rd,n))
                muc = pkd/uy
                if ( abs(muc - mu) < 1.D-2 ) exit
                uy = uy*(muc/mu)**0.5
                fy = k0*uy
                para(2) = fy
            end do
        else
            para = [ k0,pkf/mu,k1,0.D0,0.D0,0.D0,0.D0,0.D0 ]
            call rnmknl(acc,n,dt,zeta,P,ra,rv,rd,para,npara,model,tol,maxiter)
        end if

        return
        
    end subroutine rnl

    subroutine spectrum(acc,n,dt,zeta,P,np,SPA,SPI,SM) bind(c)
        integer, intent(in) :: n, np, SM
        real(8), intent(in) :: dt, zeta
        real(8), intent(in) :: acc(n), P(np)
        real(8), intent(out) :: SPA(np)
        integer, intent(out) :: SPI(np)

        select case (SM)
            case (1)
                call spafreq(acc,n,dt,zeta,P,nP,SPA,SPI)
            case (2)
                call spanmk(acc,n,dt,zeta,P,nP,SPA,SPI)
            case (3)
                call spamixed(acc,n,dt,zeta,P,nP,SPA,SPI)
            case (4)
                call pspafreq(acc,n,dt,zeta,P,nP,SPA,SPI)
            case (5)
                call pspanmk(acc,n,dt,zeta,P,nP,SPA,SPI)
            case (6)
                call pspamixed(acc,n,dt,zeta,P,nP,SPA,SPI)
            case default
                call spamixed(acc,n,dt,zeta,P,nP,SPA,SPI)
        end select

        SPA = abs(SPA)

        return

    end subroutine spectrum

    subroutine spectrumavd(acc,n,dt,zeta,P,np,SPA,SPI,SPV,SPD,SPE,SM) bind(c)
        
        integer, intent(in) :: n, np, SM
        real(8), intent(in) :: dt,zeta
        real(8), intent(in) :: acc(n), P(np)
        real(8), intent(out) :: SPA(np), SPV(np), SPD(np), SPE(np)
        integer, intent(out) :: SPI(np)

        select case (SM)
            case (1)
                call spavdfreq(acc,n,dt,zeta,P,nP,SPA,SPI,SPV,SPD,SPE)
            case (2)
                call spavdnmk(acc,n,dt,zeta,P,nP,SPA,SPI,SPV,SPD,SPE)
            case (3)
                call spavdmixed(acc,n,dt,zeta,P,nP,SPA,SPI,SPV,SPD,SPE)
            case default
                call spavdmixed(acc,n,dt,zeta,P,nP,SPA,SPI,SPV,SPD,SPE)
        end select

        SPA = abs(SPA)
        SPV = abs(SPV)
        SPD = abs(SPD)

        return

    end subroutine spectrumavd

    subroutine spmu(acc,n,dt,zeta,P,nP,SPA,SPI,SPV,SPD,SPE,mu,model,rtol,maxiter,uy,rk) bind(c)

        integer, intent(in) :: n, nP, model, maxiter
        real(8), intent(in) :: mu, dt, zeta, P(nP), rtol, rk
        real(8), intent(in) :: acc(n)
        real(8), intent(out) :: SPA(nP), SPV(nP), SPD(nP), SPE(nP)
        integer, intent(out) :: SPI(nP)
        real(8), intent(out) :: uy(nP)

        real(8) :: para(8), w, k0, k1, fy, mu1, pkd, tol
        integer :: npara, i, j, pki

        do i = 1, nP, 1
            call newmark(acc,n,dt,zeta,P(i),SPA(i),SPI(i),SPV(i),SPD(i),SPE(i))
            uy(i) = abs(SPD(i))/mu
            tol = rtol*abs(SPD(i))

            w = pi2/P(i)
            k0 = w*w
            k1 = k0*rk
            fy = k0*uy(i)
            npara = 8
            para = [ k0,fy,k1,0.D0,0.D0,0.D0,0.D0,0.D0 ]

            do j = 1, 100, 1
                call newmark_nl(acc,n,dt,zeta,P(i),SPA(i),SPI(i),SPV(i),SPD(i),SPE(i),para,npara,model,tol,maxiter)
                mu1 = abs(SPD(i))/uy(i)

                if (abs(mu1-mu)<1.D-2) exit

                uy(i) = uy(i)*(mu1/mu)**0.5
                fy = k0*uy(i)
                para(2) = fy
            end do

            SPA(i) = abs(SPA(i))
            SPV(i) = abs(SPV(i))
            SPD(i) = abs(SPD(i))
        end do

        return
    end subroutine spmu

    subroutine fitspectrum(acc,n,dt,zeta,P,nP,SPAT,a,tol,mit,fm) bind(c)

        integer, intent(in) :: n, nP, mit, fm
        real(8), intent(in) :: dt, zeta, tol
        real(8), intent(in) :: acc(n), P(nP), SPAT(nP)
        real(8), intent(out) :: a(n)

        real(8), allocatable :: EP(:), ESPAT(:)

        if ( fm == 0 ) then
            allocate(EP(nP+2))
            allocate(ESPAT(nP+2))

            EP(2:nP+1) = P
            EP(1) = P(1)*0.5
            EP(nP+2) = P(nP)*1.5
            ESPAT(2:nP+1) = SPAT
            ESPAT(1) = SPAT(1)-(SPAT(2)-SPAT(1))/(P(2)-P(1))*P(1)*0.5d0
            ESPAT(nP+2) = SPAT(nP)+(SPAT(nP)-SPAT(nP-1))/(P(nP)-P(nP-1))*P(nP)*0.5d0

            call fitspectra(acc,n,dt,zeta,EP,nP+2,ESPAT,a,tol,mit)

            deallocate(EP)
            deallocate(ESPAT)
        else if ( fm == 1 ) then
            call adjustspectra(acc,n,dt,zeta,P,nP,SPAT,a,tol,mit)
        end if
      
    end subroutine fitspectrum

end module utils